Meeting report

Alain Danet
mars 16 2022

Abstract/summary (so far)

Introduction

Why do we want to assess community changes ?

Documenting community changes

Research gap

Fish freshwater in streams

Methods

Data selection

Show code
loc <- filtered_dataset$location %>%
  left_join(filtered_dataset$site_quali, by = "siteid") %>%
  st_as_sf(coords = c("longitude", "latitude"),
  crs = 4326)
Show code
fig_sites <- paste0(
  "Australia absence because span is superior to 10 years of 100s of sites but number of samplings < 10.")
Show code
ggplot(data = world) +
    geom_sf() +
    geom_sf(data = loc, aes(color = protocol), shape = 1) +
#    coord_sf(xlim = bb[c("xmin", "xmax")],
#      ylim = bb[c("ymin", "ymax")],
#      datum = 4326, expand = TRUE) +
    theme(legend.position = "bottom") +
    labs(title = paste0("Number of sites: ", nrow(loc)))
Australia absence because span is superior to 10 years of 100s of sites but number of samplings < 10.

Figure 1: Australia absence because span is superior to 10 years of 100s of sites but number of samplings < 10.

Biodiversity

Total abundance

Total abundance was reported in four different units: Count, Count by 100m^2, Leslie index & CPUE. The two first units account for 80% of the data.

Species richness

We computed species richness with coverage-based correction proposed by (???). The Chao index correct the species richness for the sampling coverage. The sampling coverage of each sampling is evaluated by the total number of individual and the number of singleton and (double) species. All the sampling events were set at the same coverage of .80, interpolation and rarefaction for respectively low and high coverage sampling. Chao index was logged to express the temporal trends of species richness of percentage variation. The temporal trends estimated from uncorrected species richness gave the same results (Figure SXX).

Betadiversity

Betadiversity metrics compares the composition of two communities. In our case, we compare the composition of each time step relatively to the baseline, i.e. the first year of sampling.

Turnover metrics based on presence/absence can be decomposed in appearance/disappearance and nestedness / turnover.

Environment

Show code
var_pca <- dimnames(pca_riv_str$rotated$loadings)[[1]]
clean_var_pca <- get_river_atlas_significant_var()[
  get_river_atlas_significant_var() %in% var_pca]
high_load_rc1 <- pca_riv_str$rotated$loadings[
pca_riv_str$rotated$loadings[, "RC1"] > .90,
  "RC1"] 
clean_high_load_rc1_var <- get_river_atlas_significant_var()[
  get_river_atlas_significant_var() %in% names(high_load_rc1)]
Show code
fig_cap_pca <- "Rotated PCA on stream structure. The first axis is related to
the distance from the source, stralher order and discharge (m3/s). The second
axis is related to the altitude and slope degree of the site. Only the first
axis was included in the analysis"
Show code
p <- plot_rotated_pca(pca_riv_str)
Rotated PCA on stream structure. The first axis is related to
the distance from the source, stralher order and discharge (m3/s). The second
axis is related to the altitude and slope degree of the site. Only the first
axis was included in the analysis

Figure 2: Rotated PCA on stream structure. The first axis is related to the distance from the source, stralher order and discharge (m3/s). The second axis is related to the altitude and slope degree of the site. Only the first axis was included in the analysis

Show code
p$rotated
Show code
h_hft_cap <- "Distribution of difference between human footprint in 1993 and
2009."
Show code
hist(analysis_dataset$hft_ix_c9309_diff,
  xlab ="human footprint difference (1993-2009)")
Distribution of difference between human footprint in 1993 and
2009.

Figure 3: Distribution of difference between human footprint in 1993 and 2009.

Show code
s_hft_diff <- round(summary_distribution(analysis_dataset$hft_ix_c9309_diff, na.rm = T))

Statistical analysis

Species composition change

In each site and at each sampling point, we computed the dissimilarity (Jaccard, Appearance, Disappearance, Turnover & Nestedness) between each time step and the reference time, i.e. the first year of sampling. Regarding the reference time step, there are two different strategies. First is to include it, i.e. setting a dissimilarity of 0 at each beginning of a time series (Blowes et al. 2019). Another strategy is to exclude the first year from the analysis (Dornelas et al. 2014). The first option logically constrains the analysis toward increasing dissimilarity as the second can capture recovery from perturbation soon after the first samplings. I guess that the first option makes more sense more it removes to sensibility to community changes happening just after the references year.

However as the relationship between most dissimilarity indices and time is highly non linear (like a Michaelis-Menten), we logged the number of year since the first year of sampling (plus one) to linearize the relationship. AIC of the models with logged years were up to twice lower (Table SXX).

Model

Turnover

Jaccard dissimilarity, nestedness, turnover, appearance, disappearance and hillebrand were modelled in two ways, Beta and gaussian regression. Beta regressions are more appropriate for distribution bounded by [0,1] and are linearized by a logit link. Logit links have the cons that coefficients are less easy to interpret. So we run both gaussian and beta regressions and investigate the sensitivity of our results to modelling choice. Values of turnover were slightly transformed such as values of one and zero fitted beta regression requirements (\(x' = \dfrac{x \times (N - 1) + 0.5}{N}\)).

The model is defined like this:

\[ B_i = \beta_0Y_i + \beta_1Y_iR_i + \beta_2Y_iH_i + \epsilon \] where

\[ \begin{align} \\ \beta_{0}&=\beta_0 + u_{0i|j} + u_{0i} \end{align} \]

Show code
"biodiv_facet ~ 
0 + year/river_structure + year/human_footprint_diff +
(0 +  year | basin/site)"
#> [1] "biodiv_facet ~ \n0 + year/river_structure + year/human_footprint_diff +\n(0 +  year | basin/site)"

As temporal trends of biodiversity facets had a Michaelis-Menten shape, we linearized the trends by taking the log of year number added to 1. Models with logged year had a lower AIC, up to twice (Table SXX).

At each site, all the dissimilarity indices had a value of 0 the first year. Then we constrained the intercept to 0 and kept only random effect on the slope. Models comparison shows that the results are qualitatively the same (Appendix SXX).

Species richness

At the difference with dissimilarity index, we modelled species richness with the intercept and the main effect of ecological drivers. We modelled logged species richness with Chao index.

\[ B_i = \alpha + \beta_0Y_i + \beta_1R_i + \beta_2Y_iR_i + \beta_3H_i + \beta_4Y_iH_i + \epsilon \]

where \[ \begin{align} \\ \alpha&=\alpha_0 + u_{0i|j} + u_{0j} \\ \beta_{0}&=\beta_0 + u_{0i|j} + u_{0j} \end{align} \]

Total abundance

In the dataset, total abundance was described as raw count, counts by \(100m^2\), CPUE, and Leslie index representing respectively XX% of the sites. To account for this hererogeneity, we included the unit of measurement as a main affect explaining total abundance but also as an interaction with the year effect, in a similar way than previous studies (Klink et al. 2020). Similar to species richness modelling, we modelled the log of total abundances.

\[ B_i = \alpha + \beta_0Y_i + \beta_1U_i + \beta_2Y_iU_i + \beta_3R_i + \beta_4Y_iR_i + \beta_5H_i + \beta_6Y_iH_i + \epsilon \]

where \[ \begin{align} \\ \alpha&=\alpha_0 + u_{0i|j} + u_{0j} \\ \beta_{0}&=\beta_0 + u_{0i|j} + u_{0j} \end{align} \]

Implementation

All the statistical models were computed in R with the glmmTMB package. In order to facilitate the comparison of the effects, the reported slope coefficients were standardised as \(r_\delta = \beta \dfrac{\s_x}{\s_y}\), \(\beta\) being the unstandardised slope, \(s_x\) and \(s_y\) respectively the standard deviation of the explicative and response variable.

Estimation of fixed effects for each level of random effects

The fixed effect coefficients for each level of random effects were computed as sum of the random effects and its corresponding fixed effects, known as BLUPs (Best Linear Unbiased Predictions, coef.merMod in R). However, the associated standard errors and confidence intervals are not computable on BLUPs. As an example, it means that while we are able to compute the temporal trends as each site and basin, those temporal trends have no associated errors.

Show code
site_temporal_trends <- gaussian_re_self_c %>%
  select(response, random_site) %>%
  unnest(random_site) %>%
  select(response, siteid, log1_year_nb)
Show code
site_temporal_trends %>%
  filter(response %in% clust_var) %>%
  mutate(response = get_var_replacement()[response]) %>%
  ggplot(aes(x = log1_year_nb)) +
  geom_histogram() +
  facet_wrap(vars(response), scales = "free")

The estimated temporal trends of total abundance were corrected according to abundance units, with the estimated fixed effect from the model.

Cluster analysis

Show code
wide_site_tps <- site_temporal_trends  %>%
  pivot_wider(names_from = "response", values_from = "log1_year_nb")
Show code
site_tps_clust <- wide_site_tps %>%
  select(all_of(clust_var)) %>%
  as.data.frame()
rownames(site_tps_clust) <- wide_site_tps$siteid
Show code
kmeans_opt_nb_gr_sil_cap <- "Objective function value according to the
proportion of outliers data trimmed and the number of cluster(colored lines). Up
to 6 cluster, the proportion of trimmed data does not change the ranking of the clustering. We chose to stick with 6 clusters with a removal 0.05 of the outliers"
Show code
plot(clust_curv_site_fac_50)
Objective function value according to the
proportion of outliers data trimmed and the number of cluster(colored lines). Up
to 6 cluster, the proportion of trimmed data does not change the ranking of the clustering. We chose to stick with 6 clusters with a removal 0.05 of the outliers

Figure 4: Objective function value according to the proportion of outliers data trimmed and the number of cluster(colored lines). Up to 6 cluster, the proportion of trimmed data does not change the ranking of the clustering. We chose to stick with 6 clusters with a removal 0.05 of the outliers

Show code
kmeans_sites <- "Distribution of the sites in discriminant space, colored by
their cluster assignement (left). Stilhouette profile. Black points represent sites that have been trimmed
at the 5% threshold."
Show code
plot(discr_k6_fac_50)
Distribution of the sites in discriminant space, colored by
their cluster assignement (left). Stilhouette profile. Black points represent sites that have been trimmed
at the 5% threshold.

Figure 5: Distribution of the sites in discriminant space, colored by their cluster assignement (left). Stilhouette profile. Black points represent sites that have been trimmed at the 5% threshold.

Site were clustered based on the temporal trends of biodiversity facets using trimmed robust cluster method [REF]. Temporal trends were scaled but not centered before clustering. According to the optimal function, the optimal number of clusters was between 6 and 8 (Fig. ??). We used 25 starting positions and a maximum of 100 iterations. We trimmed 5% of outliers in the multidimentional space. We set cluster assignement as NA for sites that had an assignement improvement of 10% in another cluster. Clustering was computed with tclust R package.

Results

Show code
mod_comp_std_df <- gaussian_comp_std %>%
  as_tibble() %>%
  select(Parameter, starts_with("CI_"), starts_with("Coefficient")) %>%
  pivot_longer(-Parameter, names_to = "names", values_to = "values") %>%
  separate(col = names, into  = c("type", "response"), sep = "\\.") %>%
  pivot_wider(names_from = "type", values_from = "values") %>%
  filter(
    !Parameter %in% "(Intercept)",
    !str_detect(Parameter, "unitabundance"),
    response %in% clust_var 
  ) %>%
  mutate(response = get_var_replacement()[response],
    Parameter = str_replace_all(Parameter, "([a-z|0-9])\\s([a-z|0-9])", "\\1_\\2"),
    Parameter = str_replace_all(Parameter, get_model_term_replacement()) 
  )
Show code
mod_slp_std_cap <- paste0("Estimated fixed effects of the gaussian models.
  Coefficients were standardized by refiting the models with scaled response and
  explicative variables. Bars represent CI at 95% computed with Wald
  approximation.")
Show code
source("https://raw.githubusercontent.com/EmilHvitfeldt/emilfun/master/R/palette_scrapers.R")

pal <- palette_coolors("https://coolors.co/palette/ff0000-ff8700-ffd300-deff0a-a1ff0a-0aff99-0aefff-147df5-580aff-be0aff")
pal2 <- palette_coolors("https://coolors.co/palette/01befe-ffdd00-ff7d00-ff006d-adff02-8f00ff")


p <- mod_comp_std_df %>%
  ggplot(
    aes(y = Parameter, x = Coefficient, color = response, xmin = CI_low,
      xmax = CI_high)) +
  geom_vline(xintercept = 0) +
  geom_pointrange(position = position_dodge(width = .7)) +
  labs(x = "Standardized coefficients", colour = "Response variables") +
  theme(legend.position = "bottom") +
  guides(color = guide_legend(nrow = 4, byrow = TRUE))

withr::with_options(
  list(ggplot2.discrete.colour = pal),
  print(p + scale_color_discrete())
)
Estimated fixed effects of the gaussian models.
  Coefficients were standardized by refiting the models with scaled response and
  explicative variables. Bars represent CI at 95% computed with Wald
  approximation.

Figure 6: Estimated fixed effects of the gaussian models. Coefficients were standardized by refiting the models with scaled response and explicative variables. Bars represent CI at 95% computed with Wald approximation.

Show code
p_cl_dist <- site_cl_rm %>%
  select(-riv_str_rc1, -hft_c9309_scaled_no_center) %>%
  pivot_longer(-c(siteid, cl),
    names_to = "variable",
    values_to = "value") %>%
  mutate(variable = get_var_replacement()[variable]) %>%
  ggplot(aes(x = variable, y = value, color = as.factor(cl))) +
  geom_boxplot() +
  geom_hline(yintercept = 0) +
  theme(axis.text.x = element_text(angle = 20, hjust = 1)) +
  labs(x = "Response variable",
    y = "Scaled values",
    color = "Cluster")
Show code
cl_dist_cap <- "Distribution of variable values by cluster"
Show code
cl_label <- c(
  "1: No change",
  "2: High turnover",
  "3: Increase of species richness",
  "4: Decrease of species richness",
  "5: Decrease of abundance",
  "6: Increase of abundance"
)

p_cl_dist +
  scale_color_discrete(
    name = "Cluster",
    breaks = seq_len(6),
    labels = cl_label
  )
Distribution of variable values by cluster

Figure 7: Distribution of variable values by cluster

Show code
world_site_sf$cl <- world_site_sf$site %>%
    left_join(site_cl_rm[, c("siteid", "cl")], by = "siteid") %>%
    mutate(cl = relevel(as.factor(cl), ref = 1))

main_map <- ggplot(world_site_sf$world) +
  geom_sf(color = NA) +
  theme_void() +
  geom_sf(data = world_site_sf$cl %>%
    filter(!is.na(cl)),
  aes(color = cl))

zoom_map <- list(
  north_am = "USA",
  south_eu = c("FRA", "GBR", "HUN", "ESP", "BEL"),
  north_eu = c("FIN", "NOR", "SWE"),
  aus = "AUS"
)

zoom_map_coord <- map(zoom_map,
  ~world_site_sf$cl %>%
    filter(country %in% .x) %>%
    st_bbox()
)

zoom_map_coord_df <- map_dfr(zoom_map_coord,
  ~setNames(as.numeric(.x), names(zoom_map_coord[[1]])) %>%
    enframe(.),
  .id = "zoom") %>%
pivot_wider(names_from = "name", values_from = "value") %>%
mutate(label = LETTERS[seq_len(4)])

main_map_rect <- main_map +
  geom_rect(
    data = zoom_map_coord_df,
  aes(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax),
    fill = NA,
    colour = "black",
    size = 0.6) +
  geom_text(
    data = zoom_map_coord_df,
    aes(x = xmin, y = ymax, label = label),
    vjust = 1, hjust = 1, size = 15
  )

p_zoom_map <- map2(zoom_map_coord, zoom_map_coord_df$label,
  ~main_map_rect +
    coord_sf(
      xlim = .x[c("xmin", "xmax")],
      ylim = .x[c("ymin", "ymax")],
      expand = FALSE) +
    annotate("text", x = .x["xmin"], y = .x["ymax"],
      label = .y, size = 15, vjust = 1, hjust = 0) +
    theme(legend.position = "none")
)
Show code
te <- expand.grid(
  list(
    cluster = c(seq_len(6)),
    response = clust_var
  )
)

avg_cl_var <- site_cl_rm %>%
  select(-riv_str_rc1, -hft_c9309_scaled_no_center) %>%
  pivot_longer(-c(siteid, cl),
    names_to = "response",
    values_to = "value") %>%
  group_by(cl, response) %>%
  summarise(value = median(value), .groups = "drop") %>%
  rename(cluster = cl)

# Ajouter moyenne par cluster
te <- te %>%
  left_join(avg_cl_var, by = c("cluster", "response"))

heat_map_df <- te %>%
  pivot_wider(names_from = "response", values_from = "value")
Show code
heat_cl <- te %>%
  ggplot(aes(
      x = get_var_replacement()[response],
      y = str_replace_all(as.character(cluster),
        setNames(cl_label, seq_len(6))),
      fill = value)) +
  geom_tile() +
  hrbrthemes::theme_ipsum(base_family = "Helvetica") +
  coord_fixed() +
  scale_fill_gradient2(
    high = muted("red"),
    mid = "white",
    low = muted("blue"),
    midpoint = 0,
    space = "Lab") +
  labs(x = "", y = "") +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))
Show code
# Build heatmap
site_cl_rm$cl <- relevel(as.factor(site_cl_rm$cl), ref = 1)

cl_mod <- nnet::multinom(
  cl ~ riv_str_rc1 + hft_c9309_scaled_no_center,
  data = site_cl_rm)
#> # weights:  24 (15 variable)
#> initial  value 5735.422061 
#> iter  10 value 5026.238024
#> iter  20 value 4940.331090
#> final  value 4940.318756 
#> converged
Show code
ci_cl_mod <- exp(confint(cl_mod)) - 1
exp_cl_mod <- exp(coef(cl_mod))
z_cl_mod <- summary(cl_mod)$coefficients / summary(cl_mod)$standard.errors


s_cl_mod <- map_dfr(list(exp_coef = exp_cl_mod, z = z_cl_mod),
  ~.x[, c("riv_str_rc1", "hft_c9309_scaled_no_center")] %>%
    as.data.frame() %>%
    rownames_to_column(var = "cluster"), .id = "type"
) %>%
  pivot_longer(c("riv_str_rc1", "hft_c9309_scaled_no_center"),
    names_to = "variable", values_to = "value") %>%
  pivot_wider(names_from = "type", values_from = "value") %>%
  mutate(pc = exp_coef - 1)
Show code
p_cl_mod <- s_cl_mod %>%
  ggplot(aes(
      x = get_model_term_replacement()[variable],
      y = str_replace_all(as.character(cluster),
        setNames(cl_label, seq_len(6))),
      fill = pc 
      )) +
  geom_tile(aes(alpha = abs(z))) +
  hrbrthemes::theme_ipsum(base_family = "Helvetica") +
  scale_fill_continuous(labels = scales::percent_format(accuracy = 1)) +
  coord_fixed() +
  labs(x = "", y = "", alpha = "Z-score", fill = "Probability") +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))
Show code
map_cl_cap <- "(a) Distribution of the cluster biodiversity trends. NAs cluster
displays sites that could not be attributed clearly to one cluster (cf methods).
(b) Median of the scaled variables (but not centered) for each cluster. (c)
Results of the multinomial model. Color scale displays exponentiated
coefficients minus 1 while opacity displays absolute Z-scores. Color scale means that
a site has XX% more probability to be in cluster Y than in cluster 1.
"
Show code
p_map_cl <- plot_grid(
  plot_grid(main_map_rect, heat_cl, p_cl_mod, nrow = 1, labels = "auto"),
  plot_grid(plotlist = p_zoom_map, nrow = 1),
  nrow = 2
)
save_plot(here("doc/fig/p_map_cl.png"),
  p_map_cl, base_height = 4, nrow = 2, ncol = 3)
Show code
knitr::include_graphics(here("doc/fig/p_map_cl.png"))
(a) Distribution of the cluster biodiversity trends. NAs cluster
displays sites that could not be attributed clearly to one cluster (cf methods).
(b) Median of the scaled variables (but not centered) for each cluster. (c)
Results of the multinomial model. Color scale displays exponentiated
coefficients minus 1 while opacity displays absolute Z-scores. Color scale means that
a site has XX% more probability to be in cluster Y than in cluster 1.

Figure 8: (a) Distribution of the cluster biodiversity trends. NAs cluster displays sites that could not be attributed clearly to one cluster (cf methods). (b) Median of the scaled variables (but not centered) for each cluster. (c) Results of the multinomial model. Color scale displays exponentiated coefficients minus 1 while opacity displays absolute Z-scores. Color scale means that a site has XX% more probability to be in cluster Y than in cluster 1.

Show code
main_map_rect
Show code
heat_cl
Show code
p_cl_mod
Show code
p_zoom_map
#> $north_am

#> 
#> $south_eu

#> 
#> $north_eu

#> 
#> $aus

Show code
map(p_clust_prop,
  ~.x +
  scale_fill_discrete(
    name = "Cluster",
    breaks = c(0, seq_len(6)),
    labels = c("NA: undefined cluster", cl_label)
  )
)
#> $nb_cl

#> 
#> $prop_cl

Discussion

References

Blowes, Shane A., Sarah R. Supp, Laura H. Antão, Amanda Bates, Helge Bruelheide, Jonathan M. Chase, Faye Moyes, et al. 2019. “The Geography of Biodiversity Change in Marine and Terrestrial Assemblages.” Science 366 (6463): 339–45. https://doi.org/10.1126/science.aaw1620.

Dornelas, Maria, Nicholas J. Gotelli, Brian McGill, Hideyasu Shimadzu, Faye Moyes, Caya Sievers, and Anne E. Magurran. 2014. “Assemblage Time Series Reveal Biodiversity Change but Not Systematic Loss.” Science 344 (6181): 296–99. https://doi.org/10.1126/science.1248484.

Hillebrand, Helmut, Bernd Blasius, Elizabeth T. Borer, Jonathan M. Chase, John A. Downing, Britas Klemens Eriksson, Christopher T. Filstrup, et al. 2018. “Biodiversity Change Is Uncoupled from Species Richness Trends: Consequences for Conservation and Monitoring.” Journal of Applied Ecology 55 (1): 169–84. https://doi.org/https://doi.org/10.1111/1365-2664.12959.

Vellend, Mark, Lander Baeten, Isla H. Myers-Smith, Sarah C. Elmendorf, Robin Beauséjour, Carissa D. Brown, Pieter De Frenne, Kris Verheyen, and Sonja Wipf. 2013. “Global Meta-Analysis Reveals No Net Change in Local-Scale Plant Biodiversity over Time.” Proceedings of the National Academy of Sciences 110 (48): 19456–9.

(APPENDIX) Appendix

Show code
cor_site_tps_tmp_cap <-  "Spearman correlation among the temporal trends of
biodiversity facets. Temporal trends at each site computed as BLUPs predictions
(i.e. estimated from random effects). " 
Show code
cor_site_tps <- cor(
  as.matrix(select(wide_site_tps, -siteid)),
  method = "spearman")
corrplot::corrplot(cor_site_tps)
Spearman correlation among the temporal trends of
biodiversity facets. Temporal trends at each site computed as BLUPs predictions
(i.e. estimated from random effects).

Figure 9: Spearman correlation among the temporal trends of biodiversity facets. Temporal trends at each site computed as BLUPs predictions (i.e. estimated from random effects).

Analysis

Reproducibility

Reproducibility receipt

Show code
## datetime
Sys.time()
#> [1] "2022-03-16 21:05:13 CDT"
Show code

## repository
if(requireNamespace('git2r', quietly = TRUE)) {
  git2r::repository()
} else {
  c(
    system2("git", args = c("log", "--name-status", "-1"), stdout = TRUE),
    system2("git", args = c("remote", "-v"), stdout = TRUE)
  )
}
#>  [1] "commit a2801125a04b058acd1a304a6027e263e965be02"                            
#>  [2] "Author: Alain Danet <alain.danet@mnhn.fr>"                                  
#>  [3] "Date:   Wed Mar 16 09:09:19 2022 -0500"                                     
#>  [4] ""                                                                           
#>  [5] "    Try conceptual figure"                                                  
#>  [6] ""                                                                           
#>  [7] "M\t_targets.R"                                                              
#>  [8] "A\tdoc/fig/multifacet_concept.svg"                                          
#>  [9] "M\tdoc/xx-meeting-report.Rmd"                                               
#> [10] "origin\tgit@github.com:alaindanet/RivFishTimeBiodiversityFacets.git (fetch)"
#> [11] "origin\tgit@github.com:alaindanet/RivFishTimeBiodiversityFacets.git (push)"
Show code

## session info
sessionInfo()
#> R version 4.0.5 (2021-03-31)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Debian GNU/Linux 11 (bullseye)
#> 
#> Matrix products: default
#> BLAS:   /home/alain/.Renv/versions/4.0.5/lib/R/lib/libRblas.so
#> LAPACK: /home/alain/.Renv/versions/4.0.5/lib/R/lib/libRlapack.so
#> 
#> locale:
#>  [1] LC_CTYPE=fr_FR.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=fr_FR.UTF-8        LC_COLLATE=fr_FR.UTF-8    
#>  [5] LC_MONETARY=fr_FR.UTF-8    LC_MESSAGES=fr_FR.UTF-8   
#>  [7] LC_PAPER=fr_FR.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C       
#> 
#> attached base packages:
#> [1] parallel  stats     graphics  grDevices utils     datasets 
#> [7] methods   base     
#> 
#> other attached packages:
#>  [1] hrbrthemes_0.8.0        factoextra_1.0.7       
#>  [3] ade4_1.7-16             tclust_1.4-2           
#>  [5] ggeffects_1.0.1         report_0.5.0.9000      
#>  [7] see_0.6.8.1             correlation_0.8.0      
#>  [9] modelbased_0.7.1.1      effectsize_0.6.0.2     
#> [11] parameters_0.16.0       performance_0.8.0      
#> [13] bayestestR_0.11.5       datawizard_0.2.3       
#> [15] insight_0.15.0          easystats_0.4.3        
#> [17] glmmTMB_1.1.2.9000      inlatools_0.0.1.9001   
#> [19] INLA_21.11.22           sp_1.4-6               
#> [21] foreach_1.5.1           Matrix_1.3-2           
#> [23] future_1.21.0           slider_0.2.2           
#> [25] vegan_2.5-7             lattice_0.20-41        
#> [27] permute_0.9-5           codyn_2.0.5            
#> [29] janitor_2.1.0           viridis_0.5.1          
#> [31] viridisLite_0.3.0       cowplot_1.1.1          
#> [33] terra_1.5-21            rnaturalearthdata_0.1.0
#> [35] rnaturalearth_0.1.0     sf_1.0-6               
#> [37] scales_1.1.1            kableExtra_1.3.4.9000  
#> [39] here_1.0.1              lubridate_1.8.0        
#> [41] magrittr_2.0.1          forcats_0.5.1          
#> [43] stringr_1.4.0           dplyr_1.0.7            
#> [45] purrr_0.3.4             readr_2.1.1            
#> [47] tidyr_1.2.0             tibble_3.1.6           
#> [49] ggplot2_3.3.5           tidyverse_1.3.1        
#> [51] tarchetypes_0.3.2       rmarkdown_2.11         
#> [53] nvimcom_0.9-124         targets_0.8.1          
#> [55] conflicted_1.1.0       
#> 
#> loaded via a namespace (and not attached):
#>   [1] utf8_1.2.2          tidyselect_1.1.1    lme4_1.1-26        
#>   [4] grid_4.0.5          munsell_0.5.0       codetools_0.2-18   
#>   [7] units_0.8-0         distill_1.3         statmod_1.4.35     
#>  [10] withr_2.4.3         colorspace_2.0-0    highr_0.9          
#>  [13] knitr_1.36          rstudioapi_0.13     ggsignif_0.6.1     
#>  [16] Rttf2pt1_1.3.8      listenv_0.8.0       labeling_0.4.2     
#>  [19] emmeans_1.7.2       farver_2.0.3        rprojroot_2.0.2    
#>  [22] coda_0.19-4         parallelly_1.23.0   vctrs_0.3.8        
#>  [25] generics_0.1.0      TH.data_1.0-10      xfun_0.28          
#>  [28] R6_2.5.1            cachem_1.0.4        assertthat_0.2.1   
#>  [31] nnet_7.3-15         multcomp_1.4-16     gtable_0.3.0       
#>  [34] downlit_0.4.0       globals_0.14.0      processx_3.5.2     
#>  [37] sandwich_3.0-0      rlang_0.4.12        systemfonts_1.0.0  
#>  [40] splines_4.0.5       rstatix_0.7.0       extrafontdb_1.0    
#>  [43] TMB_1.7.20          broom_0.7.12        abind_1.4-5        
#>  [46] yaml_2.2.1          modelr_0.1.8        backports_1.2.1    
#>  [49] extrafont_0.17      tools_4.0.5         bookdown_0.24      
#>  [52] ellipsis_0.3.2      jquerylib_0.1.3     RColorBrewer_1.1-2 
#>  [55] Rcpp_1.0.6          progress_1.2.2      classInt_0.4-3     
#>  [58] ps_1.6.0            prettyunits_1.1.1   ggpubr_0.4.0       
#>  [61] zoo_1.8-8           haven_2.3.1         ggrepel_0.9.1      
#>  [64] cluster_2.1.1       fs_1.5.1            data.table_1.13.6  
#>  [67] openxlsx_4.2.3      warp_0.2.0          reprex_2.0.1       
#>  [70] mvtnorm_1.1-1       hms_1.1.1           evaluate_0.14      
#>  [73] xtable_1.8-4        rio_0.5.16          jpeg_0.1-8.1       
#>  [76] readxl_1.3.1        gridExtra_2.3       compiler_4.0.5     
#>  [79] KernSmooth_2.23-18  crayon_1.4.2        minqa_1.2.4        
#>  [82] htmltools_0.5.1.1   mgcv_1.8-34         tzdb_0.2.0         
#>  [85] DBI_1.1.1           corrplot_0.84       sjlabelled_1.1.7   
#>  [88] dbplyr_2.1.1        MASS_7.3-53.1       boot_1.3-27        
#>  [91] car_3.0-10          cli_3.1.0           adegraphics_1.0-15 
#>  [94] igraph_1.2.6        pkgconfig_2.0.3     foreign_0.8-81     
#>  [97] numDeriv_2016.8-1.1 xml2_1.3.2          svglite_2.0.0      
#> [100] bslib_0.2.4         webshot_0.5.2       estimability_1.3   
#> [103] rvest_1.0.2         snakecase_0.11.0    callr_3.7.0        
#> [106] digest_0.6.27       cellranger_1.1.0    gdtools_0.2.3      
#> [109] curl_4.3.2          nloptr_1.2.2.2      lifecycle_1.0.1    
#> [112] nlme_3.1-152        jsonlite_1.7.2      carData_3.0-4      
#> [115] fansi_0.5.0         pillar_1.6.4        fastmap_1.1.0      
#> [118] httr_1.4.2          survival_3.2-10     glue_1.5.1         
#> [121] zip_2.1.1           png_0.1-7           iterators_1.0.13   
#> [124] class_7.3-18        stringi_1.7.6       sass_0.3.1         
#> [127] latticeExtra_0.6-29 memoise_2.0.0       e1071_1.7-4

References